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Abstract 



We examine an extension to the theory of Gaussian wave packet dynam- 
ics in a one-dimensional potential by means of a sequence of time dependent 
displacement and squeezing transformations. Exact expressions for the quan- 
tum dynamics are found, and relationships are explored between the squeezed 
system, Gaussian wave packet dynamics, the time dependent harmonic oscil- 
lator, and wave packet dynamics in a Gauss-Hermite basis. Expressions are 
given for the matrix elements of the potential in some simple cases. Several 
examples are given, including the propagation of a non-Gaussian initial state 

in a Morse potential. 
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I. INTRODUCTION 



Wave packet dynamics has exposed interesting new phenomena in several fields. In 
femto-chemistry we are now able to time-resolve chemical processes and also observe 
effects such as the breakup and revival of wave packets 0. In atom optics wave packets 
are used to model matter waves and electron wave packets are seen in the dynamics of 
Rydberg atoms The numerical modelling of wave packet dynamics has been achieved by 
a number of methods 1^,0], but one of the earliest approaches was by Heller ^ who simply 
used the Ansatz of a time dependent Gaussian wave packet. This Gaussian approach is very 
useful, but it is usually an approximation and can be quite wrong, for example at turning 
points. Several improvements have been made: for example, the method of generalised 
Gaussian wave packets p| used complex classical trajectories for Gaussian wave packets, 
and the hybrid method used an expansion in terms of a grid of Gaussian wave packets. 

The idea of using a time dependent harmonic (i.e. Gauss-Hermite) basis, in the con- 
text of wave packet propagation, was put forward by Lee and Heller and Coalson and 
Karplus [0. The basis was chosen so that the lowest eigenstate matches the Heller Gaussian 
wave packet, but with the inclusion of a complete set of basis states the modelling can be 
performed accurately. This approach was generalised to multi-dimensional systems by Lee 
| |12| |. However, several possibilities for using a Gauss-Hermite basis exist: parameters for the 
dynamic basis were treated in a variational method by Kay and by Kucar and Meyer 
more recently the phase space picture was explored by M0ller and Henriksen |15|, and 
the use of a Gauss-Hermite basis with a variational treatment has been expanded by Billing 
ITgHTsI to examine non-adiabatic transitions and corrections to classical path equations. 

The approach that is used in this paper is similar, in principle, to the time dependent 



Gauss-Hermite basis of Refs. |jTO| and [|TT| ; in each case the basis follows the Heller Gaussian 
wave packet. However, the focus here is on the evolution operator and transformations asso- 
ciated with the time-dependent basis which becomes both "displaced" and "squeezed" . That 
is, unlike previous extensions to Gaussian wave packet dynamics, the system Hamiltonian is 
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transformed by displacement and squeezing in a way that removes all operator dependence 
which is quadratic or less. The result is an evolution equation which depends on a time de- 
pendent 'residual potential' which is based on the original (one-dimensional) potential with 
harmonic terms removed. This means that if the higher order derivatives of the potential 
are small, the system evolution will change relatively slowly, which allows a rapid numerical 
integration (of a set of ordinary differential equations). The method, which is in principle 
exact, is then similar to a time-dependent perturbation theory in a time dependent basis or 
interaction picture. Indeed, it can be developed as a perturbation theory in the higher order 
derivatives of the system potential about the classical motion. 

Because the evolution operator is found, the possibihty exists for applying this extended 
Gaussian wave packet method to initial states that are not Gaussian. The evolution operator 
also allows us to use the displacement and squeezing transformations to find explicit matrix 
elements of the residual potential using standard operator algebra. 

In section of this paper we set up the problem and perform a basis shift according to the 
classical dynamics. Section |ITT| examines Gaussian wave packet dynamics in this displaced 
basis, and in the original basis by two different approaches. The relationship between the two 
approaches is established. In section jTV] we establish the squeezing transformation necessary 
to map the time evolution of a Gaussian wave packet from its initial state. By using the 
same transformation for another change of basis we can find the equations for corrections to 
Gaussian wave packet dynamics. These equations are expressed in a Fock basis in section 0, 
where we also compare our results to the Gauss-Hermite basis. Some examples of useful 
matrix elements for potentials are given in section and in section |VI1| the results are 
applied to several different problems. 



II. DISPLACED BASIS 
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A. Scaling of the problem 



The problem we wish to describe is the one dimensional problem of a wave packet in a 
potential described by the Hamiltonian 



H 



~2 

2m 



+ U{x) . 



(1) 



The position and momentum co-ordinates have been denoted by x and p to distinguish them 
from scaled quantities that will shortly be introduced. The initial wave packet will be taken 
to be a Gaussian one. This is not essential, but it will simplify the treatment that will follow. 
The key feature is that a length scale, characteristic of the initial wave packet defines the 
width of the harmonic oscillator basis that we will use. For a Gaussian initial wave packet 



^^oix) 



exp 



{X-Xo}'' .poX 

— ^ 



(27ra2)V4 

which has a width ao and is located at xq with momentum po- 

We will now adopt a scaling of the problem such that we use the operators 
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which have the commutator [x,p] = i. The Schrodinger equation then reduces to 



where we use scaled time and energy. 



p 



+ U{x) 
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t = UJot 



(5) 



u = u/{huo) 



(6) 



The frequency ujq is determined by the width of the initial wave packet. It is the frequency 
of the harmonic oscillator for which the wave function (§) is a ground state wave function. 



In terms of the scaled quantities the initial wave function is now 



^'o(x) 
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when we use the appropriately scaled Xq and po- 



B. Local expansion of the potential 



The motion of a Gaussian wave packet in a harmonic potential is exactly solvable, even 
when the potential is time dependent, and we will use this to define the local basis for 
the wave packet. That is, the potential function U{x) will be expanded to second order 
about the the position of the wave packet. The dynamics of a Gaussian wave packet in this 
harmonic potential will be determined, and these dynamics will be used to define the basis 
for the full (non- Gaussian) wave packet dynamics. 

The wave packet (H) is located at the position xq, which will in general be time dependent; 
we then take its location to be given by Xo{t). If we expand the potential about this point 
we obtain. 



where the spatial derivatives are indicated with the primes. The terms in the expansion 
are not explicitly time dependent; they vary only with time through the position Xo{t). 
The potential function Ur{x, Xq) is the residual potential found after making the harmonic 
expansion. That is, it contains the higher order, cubic and above, terms in the expansion. 
The residual potential will play a central role in the non- Gaussian dynamics of the wave 
packet, and Eq. (|) serves as its definition. 



U{x) = U{xo) + U'{xo){x - xo) + 



U"{xo) 



(x - XoY + Ur{x,Xo) , 



(9) 
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C. Displacement of the basis 



In the following, we will make two basis changes in order to match the Gaussian part of 
the wave packet dynamics. The first basis transformation will be a displacement to remove 
the linear term in x from the potential in Eq. (^). The necessary displacement is simply 
Xo{t) in space and a momentum Po{t), such that the wave function is shifted to the origin. 
The new wave function will be 

M^,t) = D{-l3it))^ix,t) , (10) 

i.e. t) = D{f3(t))ipd{x,t), where D{P(t)) is the time dependent displacement operator 

DiPit)) = D-\-Pit)) = exp itpoit)x - txoit)p) , (11) 

with 

= (12) 

The potential U{x) in the Schrodinger equation (|^) will become transformed as 
D^^{l3{t))U {x)D{(3{t)) = U{xo + x) and we will then use the expansion in Eq. (|^). 

The requirement to remove the linear x term (and linear p term) from the potential means 
that after inserting Eq. (|TOD in the Schrodinger equation (^) we obtain the conditions: 

x'o=Po{t) (13) 
p', = -U'Mt)) (14) 

which are, of course, the classical equations of motion. With these conditions the linear 
term in x is lost and the Schrodinger equation now reads 



.dipdix,t) 



t + U{x,{t)) - if/'(xo(t))xo(t) + + Un{x, + x, a;o) 



ijd{x,t) . 

(15) 

The non-operator parts of Eq. ([T5|) are easily removed with a time dependent phase factor 



U{xo{t')) - -Xoit')U'{xo{t')) 



dt' 



-Xo{t')p',{t') + U{xoit')) 



dt' , (16) 



where for the second form we have used Eq. (|14[). Then if we define the displaced wave 
function with a phase shift as 



(17) 



we obtain the Schrodinger equation 



.di/jdp{x,t) 
dt 



p' U"Mt)) 
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+ Ur{Xo + X,Xo) 



ipdpix,t) . 



(18) 



Basis displacement has been of interest in the study of quantum state diffusion (QSD) p9[] 
where the non-linear dynamics create wave packet localisation. By using a displaced basis a 
reduction in computational effort is gained. However, in quantum state diffusion there is no 
strong motive for going to the next step of squeezing the basis because QSD localised wave 
packets all have the same size. In ordinary Schrodinger wave packet dynamics wave packets 
can change their widths enormously making basis squeezing desirable. 



III. GAUSSIAN WAVE PACKET DYNAMICS 

A. Heller's approach 

For completeness we include here an outline of standard Gaussian wave packet dynamics. 
Heller started with the Ansatz |^ 

^GVKP = exp ia{t){x - Xo{t)y + ipo{t){x - XQ{t)) + i'y{t) (19) 

in the original basis [here we use the scaled basis of Eq. (^)]. The normalisation is included 
in the time dependent complex parameter 7(t), and Heller introduced the parameter a{t) 
which characterises (the reciprocal of) the width of the Gaussian wave packet. The position 
xq and momentum po of the wave packet obey the classical equations of motion, exactly 
as in Eqs. ( P^ and (|Iip. By substituting the Gaussian wave packet into the Schrodinger 
equation @) with the truncated potential 



U{x) ~ U{xo) + U'{xo){x - xo) + ^^^(x - xof , (20) 

we can show that 

a' = -2a^ - U"{xo)/2 

^' = ia+poXo-E (21) 

where E is the classical energy pl/2 + U{xq). 

Thus the dynamics of an approximate Gaussian wave packet are completely defined by 



solving the differential equations ([T^), (|Tj), and (|2l|) . The result is approximate because 
Eq. ( ^OD is an approximation to Eq. (P). 

B. Gaussian wave packets in the displaced basis 

In order to establish some notation, and motivate the squeezing transformation in sec- 
tion this section gives an overview of Gaussian wave packet dynamics as found in the 
displaced basis of Eq. (|10|). Thus starting with Eq. (0), we again neglect the residual 
potential Ur to obtain 

.di)s{x,t) 



dt 



2 2 



i^s{x,t), (22) 



where k{t) = U"{xQ{t)) is a time dependent spring constant. This Schrodinger problem does 
have a known time dependent 'ground' state solution (see, for example, Ref. ||20|), which is 
not a stationary state, because of the time dependence in k{t). The 'ground' state solution 
can be formulated in terms of local classical trajectories. Using some of the notation of Ref. 
pO| we define a quantity e(t), through the equation 

e"it) = -k{t)e{t) , (23) 

which would make e follow the classical trajectory of a point close to the centre of the wave 
packet. For e we have the following, complex, initial conditions 
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e(0) = 1 

so that we have a time independent Wronskian with the value 

W = e'e* - ee'* = 2i . 



(24) 



(25) 



Then the ground state wave function takes the form [ pO[| 

1 



exp 



(26) 



as may be verified by substitution into Eq. (|22D. We note that the time dependent width 
of this wave packet is l/A/[2Im(e'/e)] which is found to be \e\/^/2 on using the Wronskian 
(p5|). The ground state wave function (^61) is identical to the Heller Gaussian wave packet if 
we transform it back to the original basis using the inverse of Eq. (p!7|) . That is, 



(27) 



If we perform the displacement of the wave packet we obtain Eq. with the identifications: 



a(t)=e'(t)/[2e(t)] 

7(t) = xoPo/2 - (puit) + 9n 

where Qn is a complex term arising from the normalisation of Eq. (^), i.e. 



(28) 



(29) 



Thus we see that the Heller Gaussian wave packet corresponds to the evolution of the 
'ground' state of the time dependent harmonic oscillator (E^). 



IV. TIME DEPENDENT SQUEEZED BASIS 



In order to go beyond the Gaussian wave packet approximation we need to take into 
account the non-Gaussian behaviour introduced by the residual potential Upt- This could, 
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of course, be achieved in any reasonable basis, but in order to take advantage of the power 
of Gaussian wave packet dynamics, which is often such a good approximation to the time 
evolution, it makes sense to use a time dependent basis which matches the Gaussian wave 
packet. To find this basis it is not enough to use the displaced basis of section [11 G 



we 



must also squeeze the basis to match the (time dependent) width of the Gaussian wave 
packet as well as its location. Thus, in the same way that we use a displacement operator 
to remove the linear dependence of the Hamiltonian on x in section [11 G|, we will here use 
a squeezing transformation to remove all quadratic (x and p) terms from the Hamiltonian, 
thereby leaving the naked dependence on the residual potential Ur. This will ensure that in 
the special case where the residual potential is zero, Ur = 0, the transformed wave function 
is stationary. In this case the displacement and squeezing transformations will map an 
initial Gaussian wave packet onto Heller's moving Gaussian wave packet (0). To remove 
the quadratic operator dependence we denote the squeezing transformation as IJs and define 
a wave function ipg^p in the squeezed and displaced basis as 

Writing Eq. (^ as 

. dilJdp{x,t) ^ _^ Ur{xo + X, xq)] tpdpix, t) , (31) 

we may substitute for ipdp from Eq. (|30|) to obtain 

^Ur'^^sdp + = U;' [Hs + Unixo + X, xo)] Us^sdp • (32) 

The term Hg contains all the quadratic operator dependence and can be removed if the 
operator Ug obeys 

tU7'^ = U;'HA. (33) 

To determine Ug we will start in the basis of the displaced harmonic oscillator, using the 
annihilation and creation operators 
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X + ip 
a = — 1=- 
V2 

. X — ip 
a' = — 1=— 

so that the Hamiltonian of Eq. (|2^) [Hg in Eq. (|3T|) ] becomes 

1 



2(e-e")iV-(e + e")(a' + a^'^ 



where 



= = d^d + 1/2 



(34) 



(35) 



(36) 



After some consideration we express the unitary operator Ug in the form 



Us = SiOe 



-iNd 



(37) 



where S{^) is the usual squeezing operator 



^(0 = exp L^d^' + ^-d' 



(38) 



Then using the standard expressions for the action of the squeeze operator and the phase 
shifting properties of e~*^^, the operator Us will transform the annihilation and creation 
operators as 



U;^dUs = de-'^ coshr - d^e'^'^^'^^ sinhr 
f/r^atf/, = a'fe^^coshr - de'^^^^^^ sinhr 



(39) 



where 



^ = re 



i<j> 



(40) 



Explicit expressions for r, {/), and 6, will be found by substitution in Eq. (|33D. However, to 
determine U~^^^ in Eq. ( p3| ) from the Ansatz (pTf ) we need to differentiate the exponential 
operator S with respect to the time dependence of its parameters. This is accomplished by 



first disentangling the operator, i.e. by using the relation p2 
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1 ,v 



S{^) = exp ( -^e*"^ tanh ra"^^ ) exp - ln(coshr)A^ exp ( ^'^tanhra^ ) (41) 
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1 



'«<^f oTiVi -^^2 
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(6 + ze') 






e + ie' 


(e - ^6') 



and then differentiating. It is then necessary to pull the non-exponential factors containing 
d and d^ to one side, and re-entangle the squeeze operator before a comparison can be made 
between both sides of Eq. (pSf). Some details of this calculation are presented in Appendix 
The final results for the squeezing, and phase, parameters are 

tanhr = (42) 



= -T^— 777-— 7T (43) 



= . (44) 

|e — ie\ 

A few additional relations between r, and 9 can be found in Appendix 0. 

In this way Eq. (p3D is solved and we are left with the residual potential Ur{xq + x, xq) 
in Eq. (|32D. Because Ur{xq + x, xq) depends on the operator x it will be transformed under 
the squeezing transformation (pTl). Simply by using Eq. (0) with Eq. (|39|) we find that 

f/-ix?7s = Re[e(t)]x + Im[e(t)]p . (45) 

Then we obtain from Eq. ( ^2]) the "Schrodinger equation" 

(x, t) = Ur (xo + Re[e(t)]x + Im[e(t)]p , xq) ipsdp{x, t) , (46) 



dt 

which is one of our key results. It describes the evolution of a wave function entirely in 
terms of a "potential" with higher than quadratic behaviour, i.e. in terms of the residual 
potential Ur. 

We note that neglect of Ur in Eq. (^) returns us to Gaussian wave packet dynamics. 
In this case ipsdp is stationary and we thus find from Eq. ( pO]) that 

'^Gwp{x,t) = e-''^^^''>D{P{t))Usm~\0)D~\p{0))^{x,t = 0) . (47) 

This expresses the Heller Gaussian wave packet in terms of a sequence of displacement and 
squeezing transformations, and would allow us, for example, to propagate a non-Gaussian 
wave packet in the same way as the Heller Gaussian wave packet is propagated in time. 
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V. FOCK STATE IMPLEMENTATION 



Our result so far, Eq. (^6[), describes corrections to Gaussian wave packet dynamics, but 
is hard to implement because of the appearance of the operator p throughout the transformed 
Ur. However, it is amenable to treatment in a Fock basis. If we return to the operators a 
and in Eq. (0) we can write 

U,^,U, . £W^i!M , (48) 

so that the Eq. (BUI) becomes 



z^(x,t) = Un (xo(t) + , ^^^^^^^t) . (49) 

If we now solve Eq. (|49|) in a Fock basis we obtain the non-Gaussian wave packet dynamics, 
i.e. the extended Gaussian wave packet, in the basis defined by the motion of the Gaussian 
part of the wave packet. To do this we expand the wave function in the Fock basis (of states 
labelled |n)), defined by the Gaussian wave packet ground state 

i'sdpix, t) =Y^ an(t) \n) . (50) 

n 

Then the equation of motion becomes 

m ) am{t) ■ (51) 



.danit) ^ / 



/ e*d + ea) 
\ 71 — ' ^° 



dt 

The equations which have now to be solved depend on the form of Ur and its matrix elements. 
In section ^ we will look at some specific functional forms for the residual potential in order 
to determine explicit expressions for Eq. ( ^I]) when matrix elements are taken. 

To determine the spatial wave function in the original (scaled) basis and in terms of 
the coefficients a„ we need to express Eq. ([50| ) in the original spatial basis by using the 
transformation (|30|). That is, using also Eq. (|37|), 

oo 

^(x, t) = e-^-^^W «ne-'("+'/')' {x \t){p{t))Si^) I n) . (52) 

n=0 
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We then utilize the spatial distribution of squeezed displaced Fock states and after some 
calculation obtain 



vE'(x,t) = ^^exp|z 



(X - Xo)^ + Xpo - XqPq/2 - (j)u{t) 



OO 



The values of XQ{t)^'pQ{t)^e{t\^ {t) and a^it) can be used to determine the spatial wave 
function. This result for \E' can be compared to the Ansatz employed as the starting point 
of the analysis used in Refs. , which each use a Gauss-Hermite basis. Refs. 

|l3| , p^JT7| all use a variational method where the parameters of the basis depend on the 
wave function. Refs. |lTO| , pJ| use a basis similar to Eq. (|53D, but since they chose the simplest 



kind of basis related to the Heller Gaussian wave packet ([T9| ) for n = 0, the expansion used 
differs from Eq. (|53|) by phase factors. The nth term in Eq. (|53|) has an additional phase 
of (|e|/e)". Similar phase factors, are absent in variational treatments, for example from 
Billing's Ansatz [|17| (see also Appendix ^). The n dependence of this phase factor means 
that some quadratic operator dependence is still present in the equations for the amplitudes 
of the nth Fock state. However, the variational methods try to optimise the wave packet 
trajectory — a process we do not consider here which may compensate. Also note that while 
Eq. (^) requires the evaluation of matrix elements of the residual potential, similar matrix 



elements for the full potential are required in Ref. IT^]. Again, some analytical approaches 



to these matrix elements are given in section VI 



Finally, we expect to perform a numerical integration of the various equations to deter- 
mine the wave packet dynamics of our particular system. The equations which have to be 
numerically integrated are: (i) the classical equations of motion Eqs. ([T3|), ([1^), (ii) the e 
equation (p3D, which may be split into two complex first order differential equations, or four 
real first order equations, and (iii) the amplitudes of the corrections in Eqs. (|51|). For a basis 
size of states, including the lowest energy state, this amounts to 6 + 2A^ real, first order, 
linear differential equations. The initial conditions are specified by the initial position and 
momentum of the wave packet, the initial conditions for e [Eq. (p^) ], and, in the case of a 
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Gaussian initial wave packet, oq = 1, ctn^o = 0. Then once the various matrix elements in 
Eq. ( [5l| ) have been set up, typically involving some finite sums (see the next section), the 
numerical integration can be done in a straightforward way. Note that, unlike Ref. [jl0|, it 
is not necessary to use finite difference methods on nearby trajectories, which can result in 
a reduction of numerical effort. 



VI. MATRIX ELEMENTS OF THE RESIDUAL POTENTIAL 



A. Exponential terms in Uji 



We consider a term in Ur of the form 



U/six) = exp{-(5ux) 



(54) 



where (3u is a constant characterising the potential. This term might arise from consideration 
of a Morse potential and in that case there would be two exponential terms like this one 
(see Section [VII B|) . Then, following Eq. (|49|) , we will need to evaluate 



U, 



0nin \ ^ 



' e*{t)a + e{t)a) 



m 



n 



exp 



-A 



e*{t)a + e{t)a) 



u- 



V2 



m 



(55) 



We start by disentangling the operators in the exponential 

e*(t)a + e{t)a)' 



exp 



x/2 



exp(-/3[/eaVV2) exp(-/3t/e*a/V2) exp(/?^|e|74) 



(56) 



We can then proceed several ways. For example, writing the exponentials as a power series 
and using 



a!'\m) 



{m — k)\ 



\m — k) 



(57) 



we obtain 



exp(-/3^e*a/V2)|m) = ^ 



\m — k) , 



(58) 



k=0 



k\ \l {m — k)\ 

with a similar expression for {n\ exp(— /^(yea^/ V2). Putting both of these expressions together 
we obtain the finite sums 
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2k 



• f^^k\{k + A)\{m-k)\ 



-T] 



E 

k=0 



\V 



2k 



k\{k - A)\{n - k)\ 



where 



A = n — m 



and 



n > m 
n < m 



(59) 



(60) 



(61) 



B. Power terms in Ur 



We consider a term in f/jj of the form 



(62) 



where q is an integer. Terms hke this could arise in any Taylor series expansion of a potential. 
Again, following Eq. (^), we will need to evaluate 



Uqnm ( ^ 



'e*{t)a + e{t)a^ 



m ) = { n 



e*{t)a + e{t)a) 
72 



m 



(63) 



One way to proceed is to recognise that the exponential operator in Eq. (^) can be written 



as 



exp 



-A 



e*{t)a + e{t)a) 



v 



e*{t)a + t{t)a) 
72 



(64) 



so that the coefficient of {—PuY/q^- in Eq- ( M) will lead to the matrix element (pSf ). Thus, 
after expanding the exponential in Eq. (^), and writing e in terms of its modulus and phase 
as e = lele*^" we find 

2^ 



U, 



iVnlml 



lAq^iAB^ min{m,((jr-A)/2) 



2g-A/2 



E 



qnm 



.^\q^iA9,_ min(n,{g+A)/2) 



2q+A/2 E 



A;!(A; + A)!(m-A;)![(g- A)/2-fc]! 

n > m, q — A > 

2k 



(65) 



k\{k-A)\{n-k)\[{q + A)/2-k]\ 

n < m, g + A>0 
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where for non-zero matrix elements we must have q even if A = n — m is even, or q odd 
if A = n — m is odd. We must also have g > 0. Note that these coefficients have a very 
simple dependence on |e| and 9^. This means that if the matrix elements are calculated on a 
computer, the (finite) sums do not need to be completely re-evaluated as e changes in time. 

We can write down a few examples for small g, which can also be determined by explicitly 
expanding Eq. (|63|) . For reference, some results are shown in Table | for the non-zero matrix 
elements when (n=0, 1, 2, 3...) 



C. Using a Taylor series for the potential 



We can now use the results of section |VIB| to determine the transformed residual potential 
when we expand the potential as a Taylor series about the classical position of the wave 
packet. 



U{x) = E 



fc=0 



k\ 



(66) 



with U^'^\x) as the kth. derivative of the potential. Then the residual potential is simply 
given by 



Ur{x, Xo) = 2^ □ [x - Xo) 



k=3 



k\ 



(67) 



and the transformed residual potential, as in Eq. (0), is 



Ur Xo 



e*a -f- ea 



Xq 



E 



k\ 



e*{t)a + e{t)a} 
71 



(68) 



\ V2 / fc=3 

The matrix elements of the last part of Eq. (|68|) , needed in Eq. (|5lD , can be determined by 
using the result found in Eq. (|65l), 



n 



Ur\Xq^ 7= — , Xo 



V2 



m 



^f^W(xo) 

/ , 7 1 ^ knm 



k=3 



k\ 



(69) 



For example, the Lennard- Jones potential 



TT I \ ^1 ^2 



(70) 
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may be expanded as 
leading to 



-1)^ (fc+11)! 1 „(fc + 5)!l 
I ^1 — m u ~ *^2- 



11! X 



5! 



(71) 



n 



\ — 71 — ' ^° 



m 



- / (fc + 11)! 1 _ (fc + 5)! 1 

fct^3 x^okl \ ' 11! 5! x[S 



Ufcnm • 

(72) 



VII. EXAMPLES 



A. Application to an exponential potential: Nal 

In this section we will apply the techniques developed in the previous sections to an 
example with a potential energy that varies exponentially with distance. Such a problem 
could arise in atom optics where mirrors made from evanescent waves [Q, or a periodic 



magnetisation [^], have an exponential dependence on distance from the mirror. However, 
the example we will consider here arises in a diatomic molecule where the energy of the 
molecule depends on the inter- nuclear separation (denoted by R here) . 



Wave packet dynamics in the Nal molecule have been much described [p6H 28[|, both 
theoretically and experimentally. Figure |I] shows the essential elements of the system. We 
consider two potential energy curves: one ionic and one covalent. The unexcited system is 
ionic, but at large separations of the two atoms there is an attractive Coulomb force which 
results in the level crossing near to 7 A separation. For this system we may use the potentials 
of Refs. and as, for example, used in Refs. p7| , p8| with 



V,{R) = A,exp{-(3,{R~Ro)) 



(73) 



for the covalent surface and 



V2iR) = iA2 + {B2/Rf) exp(-i?/p) - / R - e\\+ + A_)/2i?^ 

-Ci/R^ - 2e^X+\./R^ + (74) 



for the ionic surface. The centrifugal term in the dynamics is neglected. The values for the 
constants in Eqs. ( |73D and ([7^ ) are taken from (see Table [I^). 

In experiments on Nal (see Refs. |^) the ground state wave packet on the ionic surface 
is subjected to an ultra short laser pulse which places a wave packet on the covalent surface 
as indicated for t = on Fig. ^ On the covalent surface the wave packet is not in an 
equilibrium position, and so it it starts to move towards the crossing point. At the crossing 
the wave packet divides into two pieces, and the subsequent oscillations of the wave packet 
in the upper adiabatic surface have been much discussed |2^^8i. 

Here we will suppose the laser pulse is sufficiently short that its amplitude can be consid- 
ered to be a delta function. In that case the wave packet on the covalent surface at t = is 
simply proportional to the ground state wave function p7|fl| . As a result the ionic potential 
(ff^) will only serve to define the initial wave packet appearing on the covalent potential: 
we take it here to be a Gaussian with a width of 0.056 A. This narrow packet spreads quite 
rapidly, making this system of interest to an extended Gaussian wave packet method. We 
concern ourselves initially with the wave packet motion, prior to the crossing point, on the 
exponential potential (|73|). 

With the potential given by Eq. ( [75D the residual potential [Eq. (|^)] is (we will now use 
X, rather than R, for the co-ordinate to maintain the notation of sections 



Ur{x, Xq) = U (xq + x) — U (xq) — U\xq)x — U" {xq)x'^ /2 . 



(75) 



Then when we move to the squeezed basis, x becomes transformed as x — Xq + x, with 
X = U~^xUs as in Eq. (^5]). When we form the matrix elements, as in Eq. (0), we obtain 



n 



( e*a + ea) ^ 

\ ^ — 71 — ' ^° 



n 



m 



Ur [x , Xq^ 

exp(— /3i(xo(t) — -Ro)) < n| exp(— /3ix)|m > —5n^r 



~2, 



+I3i < n\x\m > —f3i < n\x \m > /2 



(76) 



Then the exponential term in Eq. ( [75| ) may be evaluated by using Eq. (pUD, and the two 
power terms may be evaluated with the q = 1,2 results from Table |l|, i.e. 
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n 



Ur (x , xo) = Aiexp(-/5i(xo(t) - i?o)) 

inm ^n,m 

+ PlUlnm-P!U2nm/2\ . (77) 



X 



With these matrix elements for the residual potential, the equations to be numerically 
integrated are: Eqs. (|TB|), (|TD, and 

The results for the Nal molecule may be seen in Fig. |I| as the two wave packets on the 
right of the figure. The full, extended wave packet is furthest to the right, and the simple 
Gaussian wave packet approximation, Eq. (|19D, is seen just slightly to the left. The wave 
packets have been derived from the amplitudes a„ by simply combining the a„ with the 
spatial representation of the harmonic oscillator eigenfunctions, as in Eq. (^), and then 
squaring the result. 

The same two wave packets are shown on a larger scale in Fig. |^. The extended Gaussian 
wave packet (XGWP) is on the right, and as well as being displaced from the Gaussian wave 
packet it does not have a Gaussian shape because its base is skewed slightly to right. These 
results have been tested against a standard numerical integration of the Schrodinger equation 
using a split step fast Fourier transform method (see, for example, Ref. for a description 
of the method). The number of basis states used for the extended Gaussian wave packet 
in Fig. II was just six. However, for this example it is found that just four basis states are 
sufficient for a reasonable approximation to the wave packet. 



B. Morse potential with Gaussian and non-Gaussian initial state 

We consider next wave packet dynamics in a Morse potential model with a potential 

UMix) =D[1- exp(-/3M(x - XM))f . (78) 

Because this potential can be expressed as the sum of two exponentials and a constant, we 
can use the results of Eqs. ([76|j77| ) (for the exponential potential) to find the matrix elements 
of the residual potential in the squeezed basis: 
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n 



Ur (x , Xo) rnj = Djexp [-2f3M{Xo{t) - Xm)] [f/2/3 nm - ^n,m] 
-2exp [-l3M{Xo{t) - Xm)] [f//3 nm - Sn,m]} 
-f/^(a;o(t))f/lnm - Ulj{xQ{t))U2nm , (79) 

where Up nm and [/2/3 nm are given by Eq. (^), f/inm and U2nm. are given in Table |, and 
U'j^j and ?7j^ are the first and second derivatives of the Morse potential ( [75] ) which will be 
evaluated at the classical position xo(t). 

The parameters for the potential (see Fig. ^ have been chosen so that the initial wave 



packet broadens considerably during the time evolution. The scaled units of Section [II A| are 
used (equivalent to m = 1,^ = 1) with an initial wave packet width of l/\/2. The results 
shown in Figures ^ and § show how the wave packet develops an asymmetry. In Fig. ^ the 
final wave packet (computed with 20 basis states) is compared to the Gaussian wave packet 
(p!9|). We see that the top of the extended Gaussian wave packet is shifted to the right in the 
Figure, and appears distinctly non-Gaussian when compared to the Gaussian wave packet. 
The 20 basis states used are sufficient to obtain convergence. 

It was briefiy mentioned at the end of Section |IV| that the squeezing and displacement 
transformations could be used to propagate a non-Gaussian wave packet. This is straight- 
forwardly done in the Fock basis of Section where it is simply a question of assigning the 
initial amplitudes in Eq. (^0|). Figure |^ shows the results of such a case where the initial 
state was chosen such that ai = 1 (with the remaining amplitudes set to zero) corresponding 
to the first excited state of a harmonic oscillator. In this case the spatial wave packet has 
the form: 

^{x) = ^^^^exp[-(a: - xo)V2] • (80) 



TT 



This initial state is propagated in the same potential shown in Fig. ^ and for the same time 
as the Gaussian initial wave packet was propagated in Fig. ^. The curve marked XGWP 
in Fig. ^ shows the result of the extended Gaussian wave packet propagation with 20 basis 
states. The dashed curve (marked UNC) shows the wave packet that results when there is 
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no coupling from the n = 1 Fock state in the dynamic basis. In this case the final wave 
packet has the same form as the initial wave packet (i.e. it is still characterised by n = 1) but 
the width, position and momentum have all changed. This means that the curve marked 
UNC amounts to the same kind of approximation to the actual final wave packet (XGWP) 
in Fig. ^ as the Gaussian wave packet in Fig. ^ is to the actual (XGWP) wave packet there. 
For Fig. I, we see that, as for the Gaussian initial wave packet in Fig. ^, it is important to 
have the coupling of the residual potential. 

The ability of the extended Gaussian wave packet method to be used for such non- 
Gaussian initial states can clearly increase the applicability of this type of method. Not 
only can the ground states of anharmonic potentials be propagated, but we could also 
propagate a thermal wave packet. In this case we would separately propagate the thermally 
populated vibrational states and then add (with thermal weightings) the final probability 
distributions. 

VIII. CONCLUSION 

In this paper we have seen a description of wave packet dynamics in terms of a time 
dependent Gaussian basis. Explicit expressions have been found for the displacement and 
squeezing parameters that describe the basis, and the displacement and squeezing trans- 
formations have been used to determine analytic expressions for matrix elements of simple 
forms of potentials. By expanding a potential as a Taylor series about the classical trajec- 
tory it is possible to use the analytic expressions for the matrix elements (in a truncated 
expansion) for almost any reasonable potential [as in Eq. (|69D]. As an example, the extended 
Gaussian wave packet method was applied to the dissociation of Nal. 

The extended Gaussian wave packet (XGWP) method is good for wave packet evolution 
where the packet remains close to a Gaussian one, and the method is especially appropriate 
if there are large changes in scale during the motion (as found in the examples treated in 
section |Vli| ). In these cases we can expect the XGWP method to be faster than a numerical 
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grid propagation method, and more accurate than a plain Gaussian wave packet method. 
Whether it is faster, or how much faster the XGWP method is, will depend on a particular 
situation. In the case of the example treated in section |V11 A|, a numerical split operator 



FFT method was found to be roughly a thousand times slower than the XGWP method. 

Finally, we should note that while the idea of the Gauss- Hermite basis has been exploited 
by a number of authors, in different ways, the emphasis in this paper has been on the 
transformations involved. Only ID results have been presented, and it is not clear if the 
method extends easily to more degrees of freedom. The method may not be so good for 
collision processes where the development of large asymmetries in the wave packet can result 
in large excitation of the squeezed basis. However, the method does seem appropriate for 
dissociative processes where there are large changes in scale, and the potential does not 
change on a length scale much smaller than the wave packet. The extended Gaussian wave 
packet method presented here can also be used to propagate non-Gaussian wave packets. 
Finally, although other Gauss-Hermite methods may have a similar numerical performance, 
it is hoped that the analytic results given here may give useful insights in the future. 
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APPENDIX A: DETERMINATION OF r, (/>, AND 9 

If we differentiate Eq. (|41|) we obtain 

1^(0 = -^a^'SiO + ^SiOa' - E'e'^^'^/^e'^^Ne^'^y^ (Al) 

where 

A = e*"^ tanh r 

E = ln(coshr). (A2) 
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Then on shifting the non-exponential N term to the right we can form 



. A' cosh r 



+i 



2B'A* + A* - A' sinh^ re 



-2i(j> 



A'e sinh r cosh r — B' 



N 



as will be required for the LHS of Eq. ( p3D . We will also need 

dt 



Now if we let 



(A3) 



(A4) 



y = coshr 
z = e**^ sinh r , 



so that 



A' = {z'y - zy')/y^ 
B' = y'/y , 



the LHS of Eq. ( p3D can be written as 
lU^ = : {zy-zy) + 



(A5) 



(A6) 



{z'*y - z*y') + iN{z'z* - y'y - 9') . (A7) 



dt 2 - - ^ ■ 2 

According to Eq. (|33|) Hg will become squeezed, and on using Eqs. (|39| ) we will find 



hi- 



2{e-e")yz + ie + e"){y' + z' 
2{e-e")yz* + ie + e"){y' + z*' 

+2N (e - e"Ky' + \z\') + (e + e"){yz + yz*) 



(A8) 



Then comparing Eq. ( |A^ ) with Eq. ( |A7| ), and inspecting the coefficient of a^^ we find that 
with the solution 
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Then comparing Eq. 
some algebra, that 



and 



It then follows that 



A = e'"^ tanh r 



e + ie' 
e — ie' 



(AlO) 



with Eq. ([A 71), and inspecting the coefficient of N we find, after 



9' 



e — le 



\e — le 



(All) 



(A12) 



coshr = -(e — ie')e 



- \t — le 



(A13) 



sinhr 



--(e + ze )e" 



(A14) 



APPENDIX B: 



To connect Eqs. (^Tj) and (p3D more closely to the work of Billing we would define 
some coefficients 



(Bl) 



and then insert Eq. ( pT|) into Eq. (|5lD , using 

^(e7e)"/2+i/4 = _,^(e7e)"/2+i/4 ^ _2zIm(a)n(e7e)"/2+V4 



to obtain, 

.<9c„(t) 



Im(a)(2ra + l)c„(t) 



+ E 



n 



( , 67t)a + 6(t)at ^ 



(B2) 



m)(e7e)("-")/2c™(t) (B3) 



as the equation of motion for the coefficients c„. 
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TABLES 
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TABLE I. 


For the potential 


x'^, values of Uqnm are shown as determined from Eq. (|65|). 


Ionic 






Covalent 


^2[eV] 
C2 [eVA6] 

A+ m 

X [A3] 
P[A] 

AE [eV] 




2760 
2.398 

11.3 
0.408 
6.431 
0.3489 
2.075 


Ai [eV] 0.813 
/3i [A-i] 4.08 
i?o [A] 2.67 



TABLE II. Parameters for the potentials (l73|) and ([74|) taken directly from Ref. 128(1 ^ 
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FIGURES 
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FIG. 1. Relevant potentials for the Nal molecule with superimposed wave packets. Expressions 



for the ionic and covalent potentials are given in Eqs. (74) and (^). The wave packets shown for 
t = and t = 0.1 ps have been scaled. The wave packet at t = has been promoted from the 
ionic potential to the covalent potential by an ultra-short pulse. The two wave packets shown at 
t = 0.1 ps arise from the Gaussian wave packet approximation (left) and the extended Gaussian 
wave packet approach (right). These two wave packets are shown more clearly in Fig. 0. 
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FIG. 2. Wave packet evolution in Nal after a time of 0.1 ps. The wave packet on the left 
(GWP) is from the Gaussian wave packet approximation, Eq. (|l9|). The wave packet on the right 
(XGWP) is the extended Gaussian wave packet from Eqs. (|^), (|l4|), (p3|) , and Eqs. (|5l|) with six 
basis states. 
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FIG. 3. Morse potential showing the initial Gaussian wave packet (left) and the final time 
evolved wave packet (right) which is shown in more detail in Fig. The initial and final wave 
packets are shown together here to show the dramatic change of width; a situation where the 
extended Gaussian wave packet approach may be appropriate. The initial wave packet is located 
at X = 19.5 (with a width of l/\/2 in scaled units) and the final wave packet is shown for t = 70.0 in 
scaled units (i.e. m = 1, ?l = 1). The Morse potential, Eq. (78), has parameters D = 10, Pm = 0.02, 
and xm = 60. 
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FIG. 4. Final wave packet, determined with 20 basis states (XGWP), for the evolution on the 
Morse potential in Fig. ^. The Gaussian wave packet (GWP) appears skewed to the left. The time 
elapsed since the initial state is 70.0 in scaled units. 
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FIG. 5. Final wave packets, determined at t = 70, for an initial n = 1 wave packet located at 
X = 19.5 in the potential shown in Fig. ^ The solid curve (XGWP) shows the extended Gaussian 
wave packet calculated in a basis of 20 states. The dashed curve (UNC) shows the final wave packet 
when there is no coupling between the basis states, and the initial wave packet only changes its 
location, width, and momentum. 
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